%to read in edi files
% compute the MT parameters
%contour / interpolate them and display as a map

%latest date 18.5.3
% 4= 0.0078,0.1875;3=0.25 ,6; 2=8:192,1=8160,256
%load c:\manoj\data\sgt2.mat;
nfiles = length(SGT2);

for i = 3:nfiles,
   		SPM=SGT2(i).SPM;
         bd4 = find(SPM.head.frq<=0.1875); 
         ln=length(bd4);
       if ln>0,
          frq=SPM.head.frq(bd4);
          for j = 1:ln,
             SPMat(j,:,:)=SPM.spectra(bd4(j)).data;
             ExEx(j) = SPMat(j,4,4);
             EyEy(j) = SPMat(j,5,5);
          end;
          Lat(i)=SPM.head.Lat;
          Long(i)=SPM.head.Log;
          C1 = pcoh(SPMat,'Ex');
          C2 = pcoh(SPMat,'Ey');
          %CSite(i)=sum(C1+C2)/(2*length(bd4));
          CSite(i) = sum((1-C1.^2)+(1-C2.^2))/(2*ln);
          %CSite(i) = sum(ExEx+EyEy)/(2*ln);
       else,
          CSite(i)=0;
          fprintf('No data for site %s\n',SPM.head.site);
       end;
       
%          SS(i,:)=sprintf('%5.3f',CSite(i));
                
          clear C1 SPMat C2 ExEx EyEy;
       end;
       
for i=3:nfiles,
S=SGT2(i).SPM.head.file;
if length(SGT2(i).SPM.head.file) ==10,
STC(i,:)=S(1:4);
elseif length(SGT2(i).SPM.head.file)==12,
STC(i,:)=S([1 2 5 6]);
end;
clear S;
end;

Long(1:2)=[];
Lat(1:2)=[];
KK=find(Lat>9&Lat<14);
D=griddata(Long(KK),Lat(KK),(CSite(KK+2)),[77:0.1:79],[10:0.1:13]');
contourf([77:0.1:79],[10:0.1:13],D);
hold on;
plot(Long(KK),Lat(KK),'k.');
%text(Long(KK),Lat(KK),STC(KK+2,:))
xlabel('Longitude(Degrees)');
ylabel('Latitude(Degrees)');
%title('Frequency Band 256 - 8160 Hz');
colorbar;